# Clear Workspace
rm(list=ls()) 
# Set number display format
options(scipen=999)
options(digits = 3)


#set the working directory accordingly
setwd("~")
###################### Libraries

library(dplyr)
library(ggplot2)
library(tidyverse)
library(data.table)
library(lubridate)

MainData=fread(".\\Data.csv", colClasses = c(ADDRESS='character'))

Data_cl=MainData %>%
        dplyr::mutate(Urban_cat=ifelse(Urban==0,"Rural", "Urban")) %>%
        dplyr::mutate(FruitVegtable=(QTY_Kg15_HH+QTY_Kg16_HH+QTY_Kg17_HH+QTY_Kg18_HH+
                        QTY_Kg19_HH+QTY_Kg20_HH)/30) %>%
        dplyr::mutate(AnimalProtein=(QTY_Kg5_HH+QTY_Kg6_HH+QTY_Kg6_HH)) %>%
        dplyr::mutate(MilkEggDairies=QTY_Kg9_HH+QTY_Kg10_HH+QTY_Kg11_HH+QTY_Kg12_HH+
                        QTY_Kg13_HH+QTY_Kg14_HH) %>%
        dplyr::mutate(Province=substr(ADDRESS,2,3))


Data_cl=Data_cl %>%
          dplyr::mutate(AnimalPortein_dum=ifelse((QTY_Kg5_HH+QTY_Kg6_HH+QTY_Kg8_HH+QTY_Kg11_HH)>0,1,0)) %>%
          dplyr::mutate(DairyBean_dum=ifelse((QTY_Kg9_HH+QTY_Kg10_HH+QTY_Kg14_HH+QTY_Kg23_HH)>0,1,0) ) %>%
          dplyr::mutate(Grains_dum=ifelse((QTY_Kg2_HH)>0,1,0)) %>%
          dplyr::mutate(Fruits_dum=ifelse((QTY_Kg15_HH+QTY_Kg16_HH)>0,1,0)) %>%
          dplyr::mutate(Vegetables_dum=ifelse((QTY_Kg17_HH+QTY_Kg18_HH+QTY_Kg19_HH+QTY_Kg20_HH)>0,1,0)) %>%
          dplyr::mutate(Variety_FoodGroups=AnimalPortein_dum*5+DairyBean_dum*5+
                          Grains_dum*5+ Fruits_dum*5+
                          Vegetables_dum*5)


ylim1 = boxplot.stats(Data_cl$QTY_Kg5_HH)$stats[c(1, 5)]

tiff(".\\Figures\\Fig1_AverageMeat_IncomeQuartile.tiff",
     width=900, height=525)

Data_cl %>%
        dplyr::filter(Year>=1370) %>%
        ggplot(aes(x=as.factor(Year), y=QTY_Kg5_HH)) + 
        geom_jitter(alpha = 0.05, color = "gray", size=.5)+
        geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
        coord_cartesian(ylim = ylim1*1.05)+
        facet_wrap(~Urban_cat, ncol=1, scale="free")+
        xlab("Year") + ylab("Average meat consumption in households")+
        theme_classic()+
        scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                                  "1373" = "1994", "1374" = "1995","1375" = "1996",
                                  "1376" = "1997", "1377" = "1998","1378" = "1999",
                                  "1379" = "2000", "1380" = "2001","1381" = "2002",
                                  "1382" = "2003", "1383" = "2004","1384" = "2005",
                                  "1385" = "2006", "1386" = "2007","1387" = "2008",
                                  "1388" = "2009", "1389" = "2010","1390" = "2011",
                                  "1391" = "2012", "1392" = "2013","1393" = "2014",
                                  "1394" = "2015", "1395" = "2016","1396" = "2017",
                                  "1397" = "2018", "1398" = "2019", "1399"="2020",
                                  "1400"="2021"))+
        theme(axis.text.x = element_text(size = 8, angle = 90),
          axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg5_HH), color="blue")


dev.off() 


tiff(".\\Figures\\Fig2_AverageMeat_IncomeQuartile.tiff",
    width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg5_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average meat consumption in households")+
  theme_classic()+scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg5_HH), color="blue")

dev.off() 


## Chicken

ylim1 = boxplot.stats(Data_cl$QTY_Kg6_HH)$stats[c(1, 5)]

tiff(".\\Figures\\FigB9_AverageChicken_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg6_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average Poultry consumption in households")+
  theme_classic()+scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg6_HH), color="blue")

dev.off() 


tiff(".\\Figures\\FigB10_AverageChicken_IncomeQuartile.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg6_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average Poultry consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg6_HH), color="blue")

dev.off() 




## Fruit vegatable

ylim1 = boxplot.stats(Data_cl$FruitVegtable)$stats[c(1, 5)]

tiff(".\\Figures\\FigB7_AverageFruitVegtable_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=FruitVegtable)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average daily fuirt and vegetable consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$FruitVegtable), color="blue")

dev.off() 


tiff(".\\Figures\\FigB8_AverageFruitVegtable_IncomeQuartile.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=FruitVegtable)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average daily fuirt and vegetable consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$FruitVegtable), color="blue")

dev.off() 



############### Seafood

ylim1 = boxplot.stats(Data_cl$QTY_Kg8_HH)$stats[c(1, 5)]

tiff(".\\Figures\\FigB11_Seafood_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg8_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average seafood consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg8_HH), color="blue")

dev.off() 


tiff(".\\Figures\\FigB12_Seafood_IncomeQuartile.tiff",
     width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg8_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average seafood consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg8_HH), color="blue")

dev.off() 



######### Milk, and dairies

ylim1 = boxplot.stats(Data_cl$MilkEggDairies)$stats[c(1, 5)]

tiff(".\\Figures\\FigB13_MilkEggDairies_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=MilkEggDairies)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average dairies and eggs consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$MilkEggDairies), color="blue")

dev.off() 


tiff(".\\Figures\\FigB14_MilkEggDairies_IncomeQuartile.tiff",
     width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=MilkEggDairies)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average  dairies and eggs consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$MilkEggDairies), color="blue")

dev.off() 


################### Flour

ylim1 = boxplot.stats(Data_cl$QTY_Kg1_HH)$stats[c(1, 5)]

tiff(".\\Figures\\FigB1_Flour_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg1_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average flour consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg1_HH), color="blue")

dev.off() 


tiff(".\\Figures\\FigB2_Flour_IncomeQuartile.tiff",
     width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg1_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average flour consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg1_HH), color="blue")

dev.off() 



################### Grain

ylim1 = boxplot.stats(Data_cl$QTY_Kg2_HH)$stats[c(1, 5)]

tiff(".\\Figures\\FigB3_Grain_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg2_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average grain consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg2_HH), color="blue")

dev.off() 


tiff(".\\Figures\\FigB4_Grain_IncomeQuartile.tiff",
     width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg2_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average grain consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg2_HH), color="blue")

dev.off() 



################### Grain

ylim1 = boxplot.stats(Data_cl$QTY_Kg2_HH)$stats[c(1, 5)]

tiff(".\\Figures\\Grain_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg2_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average grain consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg2_HH), color="blue")

dev.off() 


tiff(".\\Figures\\Grain_IncomeQuartile.tiff",
     width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg2_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average grain consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg2_HH), color="blue")

dev.off() 


################### Bread

ylim1 = boxplot.stats(Data_cl$QTY_Kg3_HH)$stats[c(1, 5)]

tiff(".\\Figures\\FigB5_Bread_UrbanRural.tiff",
     width=900, height=525)

Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg3_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~Urban_cat, ncol=1, scale="free")+
  xlab("Year") + ylab("Average bread consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg3_HH), color="blue")

dev.off() 


tiff(".\\Figures\\FigB6_Bread_IncomeQuartile.tiff",
     width=900, height=525)


Data_cl %>%
  dplyr::filter(Year>=1370) %>%
  ggplot(aes(x=as.factor(Year), y=QTY_Kg3_HH)) + 
  geom_jitter(alpha = 0.05, color = "gray", size=.5)+
  geom_boxplot(alpha = 0.1, outlier.shape=NA, color="tomato")+
  coord_cartesian(ylim = ylim1*1.05)+
  facet_wrap(~quartile, ncol=2, scale="free")+
  xlab("Year") + ylab("Average bread consumption in households")+
  theme_classic()+
  scale_x_discrete(labels=c("1370" = "1991", "1371" = "1992","1372" = "1993",
                            "1373" = "1994", "1374" = "1995","1375" = "1996",
                            "1376" = "1997", "1377" = "1998","1378" = "1999",
                            "1379" = "2000", "1380" = "2001","1381" = "2002",
                            "1382" = "2003", "1383" = "2004","1384" = "2005",
                            "1385" = "2006", "1386" = "2007","1387" = "2008",
                            "1388" = "2009", "1389" = "2010","1390" = "2011",
                            "1391" = "2012", "1392" = "2013","1393" = "2014",
                            "1394" = "2015", "1395" = "2016","1396" = "2017",
                            "1397" = "2018", "1398" = "2019", "1399"="2020",
                            "1400"="2021"))+
  theme(axis.text.x = element_text(size = 8, angle = 90),
        axis.text.y = element_text(size = 8, angle = 90))+
  annotate("rect", xmin = as.numeric(22), xmax = as.numeric(25), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  annotate("rect", xmin = as.numeric(28), xmax = as.numeric(31), ymin = 0, ymax = 15,
           alpha = .2, fill="cyan")+
  geom_hline(yintercept = mean(Data_cl$QTY_Kg3_HH), color="blue")

dev.off() 

### Generate the regression data

write.csv(Data_cl,".\\Data_cl_Regressions.csv", row.names = FALSE)

###

